import math
from decimal import Decimal
from models import ProjectsModel, SpeciesModel, SiteModel
from flask import json

def convert_decimals_to_float(obj):
    """递归地将对象中的 Decimal 转换为 float"""
    if isinstance(obj, Decimal):
        return float(obj)
    elif isinstance(obj, dict):
        return {k: convert_decimals_to_float(v) for k, v in obj.items()}
    elif isinstance(obj, list):
        return [convert_decimals_to_float(item) for item in obj]
    elif hasattr(obj, '__dict__'):  # 处理 SQLAlchemy 模型对象
        return type(obj)(**{k: convert_decimals_to_float(v)
                           for k, v in obj.__dict__.items()
                           if not k.startswith('_')})
    return obj
def expF(x: float, g0: float, gx: float, tg: float, ng: float) -> float:
    if tg != 0:
        return gx + (g0 - gx) * math.exp(-math.log(2) * (x / tg) ** ng)
    else:
        return gx
def gammaFoliage(age: float, gammaF0: float, gammaF1: float, tgammaF: float) -> float:
    if tgammaF * gammaF1 == 0:
        return gammaF1
    else:
        kgammaF = 12 * math.log(1 + gammaF1 / gammaF0) / tgammaF
        return (gammaF1 * gammaF0) / (gammaF0 + (gammaF1 - gammaF0) * math.exp(-kgammaF * age))
def lookup_table(x, x_table, y_table,interpolateLookups):
    if interpolateLookups == 'yes':
        n_table = len(x_table)
        if x <= x_table[0]:
            return y_table[0]
        elif x >= x_table[n_table - 1]:
            return y_table[n_table - 1]
        else:
            i = 0
            while x > x_table[i]:
                i += 1
            slope = (y_table[i] - y_table[i - 1]) / (x_table[i] - x_table[i - 1])
            return y_table[i - 1] + slope * (x - x_table[i - 1])
    else:
        n_table = len(x_table)
        i = 0
        while i < n_table - 1 and x > x_table[i]:
            i += 1
        return y_table[i]

def get_day_length(lat: float, day_of_year: int) -> float:
    sin_lat = math.sin(math.pi * lat / 180)
    cos_lat = math.cos(math.pi * lat / 180)
    sin_dec = 0.4 * math.sin(0.0172 * (day_of_year - 80))
    cos_h0 = -sin_dec * sin_lat / (cos_lat * math.sqrt(1 - sin_dec ** 2))

    if cos_h0 > 1:
        return 0
    elif cos_h0 < -1:
        return 1
    else:
        return math.acos(cos_h0) / math.pi


def compute_day_lengths(lat: float) -> list:
    """计算 12 个月的日长（秒数），假设一年气候数据相同"""
    days_in_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    day_of_year = -15
    day_lengths = []

    for days in days_in_month:
        day_of_year += days  # 累加计算代表日
        day_length = get_day_length(lat, day_of_year) * 86400  # 转换成秒
        day_lengths.append(day_length)

    return day_lengths


def getConductance(LAI, LAIgcx, MinCond, MaxCond, PhysMod, fCg):
    """
    计算气孔导度 (Conductance)
    """
    if LAI <= LAIgcx:
        gC = MinCond + (MaxCond - MinCond) * LAI / LAIgcx
    else:
        gC = MaxCond
    return gC * PhysMod * fCg


def getTranspiration(Q, VPD, h, gBL, gC, Qa, Qb):
    """
    计算冠层蒸腾量（kg/m²/day -> mm/day）
    """
    # 常数
    e20 = 2.2  # 20°C 时饱和水汽压随温度的变化率
    rhoAir = 1.2  # 空气密度 (kg/m³)
    lambda_ = 2460000  # 水汽化潜热 (J/kg)
    VPDconv = 0.000622  # VPD 转换系数 = 18/29/1000

    netRad = Qa + Qb * (Q * 10 ** 6 / h)  # Q 从 MJ/m²/day 转换为 W/m²

    defTerm = rhoAir * lambda_ * (VPDconv * VPD) * gBL
    div = gC * (1 + e20) + gBL
    Etransp = gC * (e20 * netRad + defTerm) / div

    transpiration = Etransp / lambda_ * h
    return transpiration


def getMortality(old_n, old_w, ms, wsx1000, thin_power):
    accuracy = 1 / 1000  # 精度，与 VBA 中的 Const accuracy 一致
    # 初始化变量
    n = old_n / 1000
    x1 = 1000 * ms * old_w / old_n
    i = 0
    while True:
        i += 1
        x2 = wsx1000 * n ** (1 - thin_power)
        fn = x2 - x1 * n - (1 - ms) * old_w
        dfn = (1 - thin_power) * x2 / n - x1
        dn = -fn / dfn
        n = n + dn

        if abs(dn) <= accuracy or i >= 5:
            break
    return old_n - 1000 * n

def to_int_month(value, default=0.0):
    MONTH_MAPPING = {
        "january": 1, "jan": 1,
        "february": 2, "feb": 2,
        "march": 3, "mar": 3,
        "april": 4, "apr": 4,
        "may": 5,
        "june": 6, "jun": 6,
        "july": 7, "jul": 7,
        "august": 8, "aug": 8,
        "september": 9, "sep": 9, "sept": 9,
        "october": 10, "oct": 10,
        "november": 11, "nov": 11,
        "december": 12, "dec": 12
    }
    if value is None:
        return default
    try:
        return int(value)
    except ValueError:
        month_str = str(value).strip().lower()
        return int(MONTH_MAPPING.get(month_str, default))
def month_name(month: int) -> str:
    if month > 12:
        month = month - 12

    if month < 0:
        month = 13

    if month == 0:
        month = 14

    match month:
        case 1:
            return "Jan"
        case 2:
            return "Feb"
        case 3:
            return "Mar"
        case 4:
            return "Apr"
        case 5:
            return "May"
        case 6:
            return "Jun"
        case 7:
            return "Jul"
        case 8:
            return "Aug"
        case 9:
            return "Sep"
        case 10:
            return "Oct"
        case 11:
            return "Nov"
        case 12:
            return "Dec"
        case 13:
            return "-"
        case 14:
            return "0"
        case _:
            return "???"

def record_monthly_result(result_list, year_month, StandAge, soil_class, FR, CO2, max_asw, min_asw,
                          DayLength, FrostDays, SolarRad, Tav, VPD, Rain, StemNo, BasArea, StandVol, avDBH,
                          MAI, MAIx, ageMAIx, ltStemGR, SLA, CanCover, LAI, LAIx, ageLAIx, WF, WR, WS, WL,
                          TotalW, AvStemMass, Density, fracBB, fT=0, fAge=0, fVPD=0, fCalpha=0, fCg=0,
                          fFrost=0, fSW=0, fNutr=0, PhysMod=0, GPP=0, NPP=0, RADint=0, alphaC=0, epsilon=0,
                          cEpsilonStem=0, CVI=0, m=0, pR=0, pS=0, pF=0, pFS=0, gammaF=0, cLitter=0, wSmax=0,
                          gammaN=0, mortality=0, SupIrrig=0, RunOff=0, fracRainIntcptn=0, RainIntcptn=0,
                          CanCond=0, WUE=0, EvapTransp=0, Transp=0, ASW=0, Height=0):
    monthly_result = {
        "YearMonth": year_month,
        "StandAge": StandAge,
        "SoilClass": soil_class,
        "FR": FR,
        "CO2": CO2,
        "maxASW": max_asw,
        "minASW": min_asw,
        "DayLength": DayLength,
        "FrostDays": FrostDays,
        "SolarRad": SolarRad,
        "Tav": Tav,
        "VPD": VPD,
        "Rain": Rain,
        "StemNo": StemNo,
        "BasArea": BasArea,
        "StandVol": StandVol,
        "avDBH": avDBH,
        "MAI": MAI,
        "MAIx": MAIx,
        "ageMAIx": ageMAIx,
        "ltStemGR": ltStemGR,
        "SLA": SLA,
        "CanCover": CanCover,
        "LAI": LAI,
        "LAIx": LAIx,
        "ageLAIx": ageLAIx,
        "WF": WF,
        "WR": WR,
        "WS": WS,
        "WL": WL,
        "TotalW": TotalW,
        "AvStemMass": AvStemMass,
        "Density": Density,
        "fracBB": fracBB,
        "fTemp": fT,
        "fAge": fAge,
        "fVPD": fVPD,
        "fCalpha": fCalpha,
        "fCg": fCg,
        "fFrost": fFrost,
        "fSW": fSW,
        "fNutr": fNutr,
        "PhysMod": PhysMod,
        "GPP": GPP,
        "NPP": NPP,
        "RadInt": RADint,
        "alphaC": alphaC,
        "Epsilon": epsilon,
        "StemEpsilon": cEpsilonStem,
        "CVI": CVI,
        "m": m,
        "pR": pR,
        "pS": pS,
        "pF": pF,
        "pFS": pFS,
        "gammaF": gammaF,
        "Litter": cLitter,
        "wSmax": wSmax,
        "gammaN": gammaN,
        "Mortality": mortality,
        "supIrrig": SupIrrig,
        "RunOff": RunOff,
        "fRainInt": fracRainIntcptn,
        "RainInt": RainIntcptn,
        "CanCond": CanCond,
        "WUE": WUE,
        "ET": EvapTransp,
        "Transp": Transp,
        "ASW": ASW,
        "Height": Height
    }
    result_list.append(monthly_result)

def calcsinglesite(projectId):
    # 获取前端传入的参数
    Pi = 3.141592654
    opStart = 0
    project = ProjectsModel.query.get(projectId)
    data = json.loads(project.data_input)
    sitename = data.get("Site", "Unknown")
    speciename = data.get("Species", "Unknown")
    latitude = data.get("Latitude", None)
    FR = data.get("Fertility_rating", None)
    soil_class_mapping = {"S": 1,"SL": 2,"CL": 3,"C": 4,"0": 0,"?": "?" }
    soil_class_str = data.get("Soil_class", "Unknown").strip().upper()
    soil_class = soil_class_mapping.get(soil_class_str, None)
    max_asw = data.get("Maximum_ASW", None)
    min_asw = data.get("Minimum_ASW", None)
    CO2 = data.get("Atmospheric_CO2", None)
    initial_wf = data.get("Initial_WF", None)
    initial_wr = data.get("Initial_WR", None)
    initial_ws = data.get("Initial_WS", None)
    stand_mass = data.get("Stand_mass", None)
    wf_fraction = data.get("WF_fraction", 0.5)
    wr_fraction = data.get("WR_fraction", 0.25)
    ws_fraction = data.get("WS_fraction", 0.25)
    if wf_fraction == 0:
        wf_fraction = 0.5
    if wr_fraction == 0:
        wr_fraction = 0.25
    if ws_fraction == 0:
        ws_fraction = 0.25
    seedling_mass = data.get("Seedling_mass", None)
    initial_date = data.get("Initial_date", None)
    initial_age_y = data.get("Initial_age_y", None)
    initial_age_m = data.get("Initial_age_m", None)
    year_planted = int(data.get("Year_planted", 0))
    month_planted = to_int_month(data.get("Month_planted", 0))
    initial_stocking = data.get("Initial_stocking", None)
    initial_asw = data.get("Initial_ASW", None)
    end_age = data.get("End_age", None)
    interpolate_lookups = data.get("Interpolate_lookups", None)

    fertility = data.get("Fertility", [])
    irrigation = data.get("Irrigation", [])
    thinning = data.get("Thinning", [])
    defoliation = data.get("Defoliation", [])
    if fertility:
        nFR = len(data["Fertility"])
        FRages = [item["Age"] for item in data["Fertility"]]
        FRVals = [item["FR"] for item in data["Fertility"]]
    else:
        nFR = 0
        FRages = []
        FRVals = []
    if irrigation:
        nIrrig = len(data["Irrigation"])
        irrigAges = [item["Age"] for item in data["Irrigation"]]
        irrigVals = [item["Ml_ha_yr"] for item in data["Irrigation"]]
    else:
        nIrrig = 0
        irrigAges = []
        irrigVals = []
    if thinning:
        nthin = len(data["Thinning"])
        thinAges = [item["Age"] for item in data["Thinning"]]
        thinVals = [item["Stocking"] for item in data["Thinning"]]
        thinWF = [item["F"] for item in data["Thinning"]]
        thinWR = [item["R"] for item in data["Thinning"]]
        thinWS = [item["S"] for item in data["Thinning"]]
    else:
        nthin = 0
        thinAges = []
        thinVals = []
    if defoliation:
        ndefol = len(data["Defoliation"])
        defolAges = [item["Age"] for item in data["Defoliation"]]
        defolVals = [item["FR"] for item in data["Defoliation"]]
    else:
        ndefol = 0
        defolAges = []
        defolVals = []
    species = convert_decimals_to_float(SpeciesModel.query.filter_by(speciesname=speciename).first())
    site = convert_decimals_to_float(SiteModel.query.filter_by(sitename=sitename).first())
    # 提取站点气候数据
    climate_data = {
        "tmax": [],
        "tmin": [],
        "precip": [],
        "evap": [],
        "ra": [],
        "raind": [],
        "frost": []
    }
    for month in range(1, 13):
        tmax = getattr(site, f'tmax{month}', None)
        tmin = getattr(site, f'tmin{month}', None)
        precip = getattr(site, f'precip{month}', None)
        evap = getattr(site, f'evap{month}', None)
        ra = getattr(site, f'ra{month}', None)
        raind = getattr(site, f'raind{month}', None)
        frost = getattr(site, f'frost{month}', None)
        climate_data["tmax"].append(tmax)
        climate_data["tmin"].append(tmin)
        climate_data["precip"].append(precip)
        climate_data["evap"].append(evap)
        climate_data["ra"].append(ra)
        climate_data["raind"].append(raind)
        climate_data["frost"].append(frost)
    climate_data["daylength"] = compute_day_lengths(latitude)
    poolFractn = 0
    pFS2 = species.pFS2
    pFS20 = species.pFS20
    aWs = species.aS
    nWs = species.nS
    pRx = species.pRx
    pRn = species.pRn
    gammaFx = species.gammaFx
    gammaF0 = species.gammaF0
    tgammaF = species.tgammaF
    gammaR = species.gammaR
    Tmin = species.Tmin
    Topt = species.Topt
    Tmax = species.Tmax
    kF = species.kF
    fCalpha700 = species.fCalpha700
    fCg700 = species.fCg700
    m0 = species.m0
    fN0 = species.fN0
    fNn = species.fNn
    MaxAge = species.MaxAge
    nAge = species.nAge
    rAge = species.rAge
    gammaNx = species.gammaNx
    gammaN0 = species.gammaN0
    tgammaN = species.tgammaN
    ngammaN = species.ngammaN
    wSx1000 = species.wSx1000
    thinPower = species.thinPower
    mF = species.mF
    mR = species.mR
    mS = species.mS
    SLA0 = species.SLA0
    SLA1 = species.SLA1
    tSLA = species.tSLA
    k = species.k
    fullCanAge = species.fullCanAge
    MaxIntcptn = species.MaxIntcptn
    LAImaxIntcptn = species.LAImaxIntcptn
    alphaCx = species.alpha
    Y = species.Y
    MinCond = species.MinCond
    MaxCond = species.MaxCond
    LAIgcx = species.LAIgcx
    CoeffCond = species.CoeffCond
    BLcond = species.BLcond
    fracBB0 = species.fracBB0
    fracBB1 = species.fracBB1
    tBB = species.tBB
    rho0 = species.rhoMin
    rho1 = species.rhoMax
    tRho = species.tRho
    aH = species.aH
    nHB = species.nHB
    nHN = species.nHN
    aV = species.aV
    nVB = species.nVB
    nVN = species.nVN
    Qa = species.Qa
    Qb = species.Qb
    gDM_mol = species.gDM_mol
    molPAR_MJ = species.molPAR_MJ
    # print("climate_data:", climate_data)
    result = []

    if soil_class > 0:
        SWconst = 0.8 - 0.1 * soil_class
        SWpower = 11 - 2 * soil_class
    elif soil_class == 0:
        SWconst = 999
        SWpower = species.SWpower
    else:
        SWconst = species.SWconst
        SWpower = species.SWpower
    if fNn == 0:
        fN0 = 1
    pfsPower = math.log(pFS20 / pFS2) / math.log(20 / 2)
    pfsConst = pFS2 / (2 ** pfsPower)
    fCalphax = fCalpha700 / (2 - fCalpha700)
    fCg0 = fCg700 / (2 * fCg700 - 1)
    if min_asw > max_asw:
        min_asw = max_asw
    if initial_asw <= min_asw:
        initial_asw = min_asw
    if initial_asw >= max_asw:
        initial_asw = max_asw
    poolFractn = max(0, min(1, poolFractn))
    applIrrig = 0
    pooledSW = 0
    MAIx = 0
    LAIx = 0
    daysInMonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    interpolateLookups = 'no' #options里面的!!!!!!!!
    opFrequency = 3 # 输出频率！！！！！！！！
    if opFrequency == 3:
        monthsInStep = 1
    else:
        monthsInStep = 12
    opDateNo = 1
    if seedling_mass>0:
        initial_year = year_planted
        initial_month = month_planted
    elif initial_age_y!=0 or initial_age_m!=0:
        initial_year = year_planted + initial_age_y
        initial_month = month_planted + initial_age_m
        while initial_month>=12:
            initial_year += 1
            initial_month -= 12
    else:
        initial_year = int(data.get("Initial_year", 0))
        initial_month = to_int_month(data.get("Initial_month", 0))
        if initial_year < year_planted:
            initial_year = year_planted + initial_year
    if seedling_mass > 0:
        initial_wr = (wr_fraction / (wf_fraction + wr_fraction + ws_fraction)) * seedling_mass * initial_stocking / (10**6)
        initial_wf = (wf_fraction / (wf_fraction + wr_fraction + ws_fraction)) * seedling_mass * initial_stocking / (10**6)
        initial_ws = (ws_fraction / (wf_fraction + wr_fraction + ws_fraction)) * seedling_mass * initial_stocking / (10**6)
    elif stand_mass > 0:
        initial_wr = (wr_fraction / (wf_fraction + wr_fraction + ws_fraction)) * stand_mass
        initial_wf = (wf_fraction / (wf_fraction + wr_fraction + ws_fraction)) * stand_mass
        initial_ws = (ws_fraction / (wf_fraction + wr_fraction + ws_fraction)) * stand_mass
    StandAge = (initial_year + initial_month / 12) - (year_planted + month_planted / 12)
    StartAge = int(StandAge)
    WS = initial_ws
    WF = initial_wf
    WR = initial_wr
    WL = 0
    TotalW = initial_ws + initial_wf + initial_wr
    StemNo = initial_stocking
    ASW = initial_asw
    thinEventNo = 1
    defoltnEventNo = 1
    ltStemGR = 0
    ltStartAge = 4
    if StandAge > ltStartAge:
        ltAge0 = StandAge
    else:
        ltAge0 = ltStartAge
    SLA = expF(StandAge, SLA0, SLA1, tSLA, 2)
    fracBB = expF(StandAge, fracBB0, fracBB1, tBB, 1)
    Density = expF(StandAge, rho0, rho1, tRho, 1)
    gammaF = gammaFoliage(StandAge,gammaF0, gammaFx, tgammaF)

    AvStemMass = WS * 1000 / StemNo
    avDBH = (AvStemMass / aWs) ** (1 / nWs)
    BasArea = (((avDBH / 200) ** 2) * Pi) * StemNo
    Height = aH * avDBH ** nHB * StemNo ** nHN
    LAI = WF * SLA * 0.1
    cLAI = LAI

    if aV > 0:
        StandVol = aV * (avDBH ** nVB) * (StemNo ** nVN)
    else:
        StandVol = WS * (1 - fracBB) / Density
    oldVol = StandVol
    if StandAge > 0:
        MAI = StandVol / StandAge
    else:
        MAI = 0
    initial_year_display = int(year_planted + month_planted / 12 + StandAge - 0.001)
    initial_month_name = month_name(initial_month)
    initial_year_month = f"{initial_year_display} - {initial_month_name}"
    initial_index = initial_month - 1
    record_monthly_result(result_list=result,year_month=initial_year_month,StandAge=StandAge,soil_class=soil_class,FR=FR,CO2=CO2,max_asw=max_asw,min_asw=min_asw,
        DayLength=climate_data["daylength"][initial_index],FrostDays=climate_data["frost"][initial_index],SolarRad=climate_data["ra"][initial_index],Tav=(climate_data["tmax"][initial_index] + climate_data["tmin"][initial_index]) / 2,
        VPD=((6.1078 * math.exp(17.269 * climate_data["tmax"][initial_index] / (237.3 + climate_data["tmax"][initial_index]))) - (6.1078 * math.exp(
                 17.269 * climate_data["tmin"][initial_index] / (237.3 + climate_data["tmin"][initial_index])))) / 2,Rain=climate_data["precip"][initial_index],StemNo=StemNo,
        BasArea=BasArea,StandVol=StandVol,avDBH=avDBH,MAI=MAI,MAIx=MAIx,ageMAIx=0,ltStemGR=ltStemGR,SLA=SLA,
        CanCover=1 if StandAge >= fullCanAge else (StandAge + 0.01) / fullCanAge,LAI=LAI,LAIx=LAIx,ageLAIx=0,WF=WF,WR=WR,WS=WS,WL=WL,TotalW=TotalW,AvStemMass=AvStemMass,
        Density=Density,fracBB=fracBB,gammaF=gammaF,wSmax=wSx1000 * (1000 / StemNo) ** thinPower,ASW=ASW,Height=Height
    )
    metMonth = initial_month
    for year in range(end_age - StartAge):
        # 初始化年度累计变量
        aStemDM = aRADint = aGPP = aNPP = aEvapTransp = aTransp = aRunOff = aSupIrrig = ageMAIx = ageLAIx = 0

        month = initial_month-1
        for monthCounter in range(12):
            # 初始化月度累积变量
            if opFrequency == 3 or monthCounter == 0:
                selfThin = mortality = cRADint = cLAI = cCVI = cGPP = cNPP = cStemDM = 0
                cTransp = cEvapTransp = cRainInt = cRunOff = cSupIrrig = cLitter = 0
            if abs(StandAge - ltAge0) < 0.001:
                ltWS0 = WS
            if nFR > 0:
                FR = lookup_table(StandAge, FRages, FRVals,interpolateLookups)
            if nIrrig > 0:
                applIrrig = lookup_table(StandAge, irrigAges, irrigVals,interpolateLookups)
            # 读取气象数据
            month = (month + 1) % 12
            metMonth = metMonth + 1
            if metMonth > 12:
                metMonth = 1
            index = (metMonth - 1) % 12
            SolarRad = climate_data["ra"][index]
            Tx = climate_data["tmax"][index]
            Tn = climate_data["tmin"][index]
            Tav = (Tx + Tn) / 2
            VPD =( (6.1078 * math.exp(17.269 * Tx / (237.3 + Tx))) - (6.1078 * math.exp(17.269 * Tn / (237.3 + Tn)))) / 2
            FrostDays = climate_data["frost"][index]
            Rain = climate_data["precip"][index]
            DayLength = climate_data["daylength"][index]

            # 计算环境修正因子
            fT = 0 if Tav <= Tmin or Tav >= Tmax else (((Tav - Tmin) / (Topt - Tmin)) * \
                                                       ((Tmax - Tav) / (Tmax - Topt)) ** (
                                                                   (Tmax - Topt) / (Topt - Tmin)))
            fVPD = math.exp(-CoeffCond * VPD)
            MoistRatio = ASW / max_asw
            fSW = 1 / (1 + ((1 - MoistRatio) / SWconst) ** SWpower)

            fNutr = 1 if fNn == 0 else 1 - (1 - fN0) * (1 - FR) ** fNn
            fFrost = 1 - kF * (FrostDays / 30)
            fCalpha = fCalphax * CO2 / (350 * (fCalphax - 1) + CO2)
            fCg = fCg0 / (1 + (fCg0 - 1) * CO2 / 350)

            if nAge == 0:
                fAge = 1
            else:
                fAge = 1 / (1 + (StandAge / MaxAge / rAge) ** nAge)
                RelAge = StandAge / MaxAge
            PhysMod = min(fVPD, fSW) * fAge
            if fullCanAge > 0 and StandAge < fullCanAge:
                CanCover = (StandAge + 0.01) / fullCanAge
            else:
                CanCover = 1
            lightIntcptn = 1 - math.exp(-k * LAI / CanCover)

            # 计算光合速率
            alphaC = alphaCx * fNutr * fT * fFrost * fCalpha * PhysMod
            epsilon = gDM_mol * molPAR_MJ * alphaC
            RAD = SolarRad * daysInMonth[month]
            RADint = RAD * lightIntcptn * CanCover
            GPP = epsilon * RADint / 100
            NPP = GPP * Y

            # 计算水分平衡
            CanCond = getConductance(LAI,LAIgcx, MinCond, MaxCond, PhysMod, fCg)
            Transp = daysInMonth[month] * getTranspiration(SolarRad, VPD, DayLength, BLcond, CanCond, Qa, Qb)
            if LAImaxIntcptn > 0:
                fracRainIntcptn = MaxIntcptn * min(1, LAI / LAImaxIntcptn)
            else:
                fracRainIntcptn = MaxIntcptn

            RainIntcptn = Rain * fracRainIntcptn
            SupIrrig = 0
            RunOff = 0
            ASW = ASW + Rain + (100 * applIrrig / 12) + pooledSW
            EvapTransp = min(ASW, Transp + RainIntcptn)
            excessSW = max(ASW - EvapTransp - max_asw, 0)
            ASW = ASW - EvapTransp - excessSW
            pooledSW = poolFractn * excessSW
            RunOff = (1 - poolFractn) * excessSW
            if ASW < min_asw:
                SupIrrig = min_asw - ASW
                ASW = min_asw
            TranspScaleFactor = EvapTransp / (Transp + RainIntcptn)
            GPP = TranspScaleFactor * GPP
            NPP = TranspScaleFactor * NPP
            if EvapTransp != 0:
                WUE = 100 * NPP / EvapTransp
            # 计算生物量增量
            m = m0 + (1 - m0) * FR
            pFS = pfsConst * (avDBH ** pfsPower)
            pR = pRx * pRn / (pRn + (pRx - pRn) * PhysMod * m)
            pS = (1 - pR) / (1 + pFS)
            pF = 1 - pR - pS

            incrWF = NPP * pF
            incrWR = NPP * pR
            incrWS = NPP * pS
            lossWF = gammaF * WF
            lossWR = gammaR * WR
            WF = WF + incrWF - lossWF
            WR = WR + incrWR - lossWR
            WS = WS + incrWS
            WL = WL + lossWF
            TotalW = WF + WR + WS
            StandAge += 1 / 12
            if thinEventNo <= nthin:
                if StandAge >= thinAges[thinEventNo-1]:
                    if StemNo > thinVals[thinEventNo-1]:
                        delN = (StemNo - thinVals[thinEventNo-1]) / StemNo

                        StemNo = StemNo * (1 - delN)
                        WF = WF * (1 - delN * thinWF[thinEventNo-1])
                        WR = WR * (1 - delN * thinWR[thinEventNo-1])
                        WS = WS * (1 - delN * thinWS[thinEventNo-1])
                    thinEventNo = thinEventNo + 1
            if defoltnEventNo <= ndefol:
                if StandAge >= defolAges[defoltnEventNo-1]:
                    WF = WF * defolVals[defoltnEventNo-1]
                defoltnEventNo = defoltnEventNo + 1
            gammaN = expF(StandAge, gammaN0, gammaNx, tgammaN, ngammaN)
            if gammaN > 0:
                delStems = gammaN * StemNo / 12 / 100
                WF = WF - mF * delStems * (WF / StemNo)
                WR = WR - mR * delStems * (WR / StemNo)
                WS = WS - mS * delStems * (WS / StemNo)
                StemNo = StemNo - delStems
                mortality = mortality + delStems
            wSmax = wSx1000 * (1000 / StemNo) ** thinPower
            AvStemMass = WS * 1000 / StemNo
            if wSmax < AvStemMass:
                delStems = getMortality(StemNo, WS, mS, wSx1000, thinPower)
                WF = WF - mF * delStems * (WF / StemNo)
                WR = WR - mR * delStems * (WR / StemNo)
                WS = WS - mS * delStems * (WS / StemNo)
                StemNo = StemNo - delStems
                wSmax = wSx1000 * (1000 / StemNo) ** thinPower
                AvStemMass = WS * 1000 / StemNo
                selfThin = selfThin + delStems
            SLA = expF(StandAge, SLA0, SLA1, tSLA, 2)
            fracBB = expF(StandAge, fracBB0, fracBB1, tBB, 1)
            Density = expF(StandAge, rho0, rho1, tRho, 1)
            gammaF = gammaFoliage(StandAge,gammaF0, gammaFx, tgammaF)
            LAI = WF * SLA * 0.1
            avDBH = (AvStemMass / aWs) ** (1 / nWs)
            BasArea = (((avDBH / 200) ** 2) * Pi) * StemNo
            Height = aH * (avDBH ** nHB) * (StemNo ** nHN)
            if aV > 0:
                StandVol = aV * (avDBH ** nVB) * (StemNo ** nVN)
            else:
                StandVol = WS * (1 - fracBB) / Density
            CVI = StandVol - oldVol
            oldVol = StandVol
            if StandAge > 0:
                MAI = StandVol / StandAge
            else:
                MAI = 0
            cRADint = cRADint + RADint
            aRADint = aRADint + RADint
            cGPP = cGPP + GPP
            aGPP = aGPP + GPP
            cNPP = cNPP + NPP
            aNPP = aNPP + NPP
            cCVI = cCVI + CVI
            cLitter = cLitter + lossWF
            cStemDM = cStemDM + incrWS
            aStemDM = aStemDM + incrWS
            cRainInt = cRainInt + RainIntcptn
            cTransp = cTransp + Transp
            aTransp = aTransp + Transp
            cEvapTransp = cEvapTransp + EvapTransp
            aEvapTransp = aEvapTransp + EvapTransp
            cRunOff = cRunOff + RunOff
            aRunOff = aRunOff + RunOff
            aSupIrrig = aSupIrrig + SupIrrig
            cSupIrrig = cSupIrrig + SupIrrig
            cLAI = cLAI + LAI / monthsInStep
            if cRADint != 0:
                cEpsilonStem = 100 * cStemDM / cRADint
            else:
                cEpsilonGross = 100 * cGPP / cRADint
            if cEvapTransp != 0:
                cWUE = 100 * cNPP / cEvapTransp
            if StandAge > ltAge0 + 0.99:
                ltStemGR = (WS - ltWS0) / (StandAge - ltAge0)
            if (opFrequency == 3) or (monthCounter == 12):
                if cLAI > LAIx:
                    LAIx = cLAI
                    ageLAIx = StandAge
                if MAI > MAIx:
                    MAIx = MAI
                    ageMAIx = StandAge
            # 写入月度结果
            current_month = (initial_month + monthCounter) % 12
            if current_month == 0:
                current_month = 12
            current_month = (initial_month + monthCounter) % 12
            if current_month == 0:
                current_month = 12
            if opFrequency == 3 or (opFrequency == 2 and current_month == initial_month):
                year = int(year_planted + month_planted / 12 + StandAge - 0.001)
                monthname = month_name(initial_month + monthCounter)
                year_month = f"{year} - {monthname}"
                record_monthly_result(
                    result_list=result,
                    year_month=year_month,
                    StandAge=StandAge,
                    soil_class=soil_class,
                    FR=FR,
                    CO2=CO2,
                    max_asw=max_asw,
                    min_asw=min_asw,
                    DayLength=DayLength,
                    FrostDays=FrostDays,
                    SolarRad=SolarRad,
                    Tav=Tav,
                    VPD=VPD,
                    Rain=Rain,
                    StemNo=StemNo,
                    BasArea=BasArea,
                    StandVol=StandVol,
                    avDBH=avDBH,
                    MAI=MAI,
                    MAIx=MAIx,
                    ageMAIx=ageMAIx,
                    ltStemGR=ltStemGR,
                    SLA=SLA,
                    CanCover=CanCover,
                    LAI=LAI,
                    LAIx=LAIx,
                    ageLAIx=ageLAIx,
                    WF=WF,
                    WR=WR,
                    WS=WS,
                    WL=WL,
                    TotalW=TotalW,
                    AvStemMass=AvStemMass,
                    Density=Density,
                    fracBB=fracBB,
                    fT=fT,
                    fAge=fAge,
                    fVPD=fVPD,
                    fCalpha=fCalpha,
                    fCg=fCg,
                    fFrost=fFrost,
                    fSW=fSW,
                    fNutr=fNutr,
                    PhysMod=PhysMod,
                    GPP=GPP,
                    NPP=NPP,
                    RADint=RADint,
                    alphaC=alphaC,
                    epsilon=epsilon,
                    cEpsilonStem=cEpsilonStem,
                    CVI=CVI,
                    m=m,
                    pR=pR,
                    pS=pS,
                    pF=pF,
                    pFS=pFS,
                    gammaF=gammaF,
                    cLitter=cLitter,
                    wSmax=wSmax,
                    gammaN=gammaN,
                    mortality=mortality,
                    SupIrrig=SupIrrig,
                    RunOff=RunOff,
                    fracRainIntcptn=fracRainIntcptn,
                    RainIntcptn=RainIntcptn,
                    CanCond=CanCond,
                    WUE=WUE,
                    EvapTransp=EvapTransp,
                    Transp=Transp,
                    ASW=ASW,
                    Height=Height
                )
        if aRADint != 0:
            aEpsilonStem = 100 * aStemDM / aRADint
            aEpsilonGross = 100 * aGPP / aRADint
        if aEvapTransp != 0:
            aWUE = 100 * aNPP / aEvapTransp

        if aRADint != 0:
            aEpsilonStem = 100 * aStemDM / aRADint
            aEpsilonGross = 100 * aGPP / aRADint
        if aEvapTransp != 0:
            aWUE = 100 * aNPP / aEvapTransp

        # 写入年度结果
        # write3PGresults(opEndYear, year, 12)

    return result